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Abstract 

We propose a new method to solve the Hartree-Fock-Bogoliubov equations for 
weakly bound nuclei whose purpose is to improve the treatment of the contin- 
uum when a finite range two-body interaction is used. We replace the traditional 
expansion on a discrete harmonic oscillator basis by a mixed eigenfunction expan- 
sion associated with a potential that explicitely includes a continuous spectrum. We 
overcome the problem of continuous spectrum discretization by using a resonance 
expansion. 
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1 Introduction 



In many problems of present-day nuclear physics (structure of the ground state 
of nuclei close to drip lines, description of weakly bound or unbound excited 
states, decay of isomeric states), one has to explicitly take into account the 
continuum of scattering states in addition to the discrete spectrum of bound 
states. A convenient approach to these problems is provided by the mean- 
field theory in which pairing correlations are included in a self-consistent way 
through the Hartree-Fock-Bogoliubov (HFB) formalism [1,2]. In this method a 
set of quasi-particle excitations together with their vacuum is obtained, which 
forms the basis for a description of ground state properties and spectra. In 
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principle, quasi-particle states, when expressed in terms of the single-particle 
eigenstates of a one-body hamiltonian, are always superpositions of states 
belonging to both the discrete and continuous spectra [3,4]. In ordinary sta- 
ble nuclei, the continuum part of quasi-particle states is usually neglected, 
as ground state properties and low-energy excitations mainly depend on the 
discrete quasi-particle states representing occupied and negative energy un- 
occupied single-particle orbitals. So, quasi-particles are often obtained using 
expansions on a discrete basis of orthonormal square-integrable functions. 

A very different situation is encountered in weakly bound nuclei, where the gap 
between the last occupied single-particle level and the continuum of positive 
energy states may become smaller than a few MeV. In this case, the residual 
pairing interaction is able to induce a significant coupling between the discrete 
and continuum single-particle states and the contribution of the continuum to 
quasi-particle states cannot be neglected. 

In such a situation, new numerical methods for solving the HFB equations 
have to be implemented, which more or less implicitly, amount to introduce 
a suitable discretization of the continuum. With these methods, not only the 
contribution of the continuum to quasi-particle states can be obtained, but 
also a description of quasi-particles belonging to the continuous spectrum, 
including resonances. Let us note that such a treatment is necessary to repro- 
duce the unusual properties of weakly bound nuclei such as halos and neutron 
skins. 

Several techniques have been proposed in the past along these lines: lattice 
calculations [5], Basis-Spline Galerkin lattice [6], methods using the local-scale 
point transformation of the spherical harmonic oscillator wave functions [7], 
which are adapted to calculations involving zero-range two-body interactions 
such as the Skyrme forces. Another approach, based on the Kamimura-Gauss 
basis [8] has been reported recently [9]. It appears to be a promising tech- 
nique for both zero range and finite range forces. Recently, a method has been 
proposed where the HFB equations are solved in a spherical box with exact 
boundary conditions for scattering states [10,11]. In these methods, the HFB 
equations are usually expressed in space coordinates, which is especially con- 
venient when a zero-range nucleon-nucleon interaction is used, and interesting 
results have been obtained in this context [3,12,13]. 

In the case of the finite range Gogny interaction [2], the HFB equations have 
always been solved by expanding the quasi-particle states on a finite discrete 
basis of orthonormal functions. In view of computational convenience, the lat- 
ter are usually chosen as the eigenfunctions of an harmonic oscillator (HO), 
whose geometry and symmetries are adapted to the nuclear states to be stud- 
ied. The HFB equations are then expressed in matrix form, and they are solved 
by using an iterative procedure. Namely, at iteration n, the set of individual 
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wave functions } is expanded on a finite basis of generalized eigenstates 
of a reference hamiltonian H (n \ generally taken as T + vjj% where is an 
harmonic oscillator potential: 

Nmax 

4 n) = E c£V?>- (i) 
P =i 



The are explicitly known and the unknown coefficients Cj™ are deter- 
mined using matrix diagonalization techniques. 

This simple choice has been found to be appropriate for computing the prop- 
erties of well bound nuclei, the convergence of the iterative procedure being 
obtained with a reasonable number of iterations [2,14]. In practice, the sum in 
(1) is truncated to a number of states N max corresponding to a number of HO 
shells ranging from 6-8 in light nuclei up to 16-20 in heavy and superheavy 
nuclei. 

On the other hand, in exotic nuclei close to drip-lines, weakly bound single- 
particle states acquire a spatial extension so large that the use of expansions 
on HO wave functions such as (1) requires prohibitively large values of N max . 
In this mentioned above, an alternative method of solving the HFB 

equations has to be found. 

A natural way to avoid this difficulty is to change the reference potential 
from which the basis is built. Let us recall the well-known "decomposition 
of unity" on a complete basis of generalized eigenstates, associated with a 
reference hamiltonian, a sum of kinetic and potential terms H = T + V. We 
assume that V is short-ranged and negative for small r. This decomposition 
formula can be written formally as 

N 

1 = E l^n >< Vn\ + / \<f(; k) >< (f(; k)\ dk, (2) 
n=l J 

where the first (resp. second) contribution corresponds to the discrete (resp. 
continuous) spectrum of H. 

Numerical simulations [15,16] show that the continuous contribution in (2) can 
be neglected in the expansion of single particle states provided the gap between 
the energy of the last occupied state \E^\ and the continuous threshold is large 
enough, which is the case in stable nuclei. 

As mentioned before, for nuclei near instability, the gap is small and the 
continuous contribution in (2) has to be taken into account. 

A new difficulty then appears: a numerical treatment based on equation (2) 
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requires to introduce both a discretization and a cut-off in the integral in- 
volving continuous states. This leads to a finite expansion upon discrete basis 
states which is not necessarly less expensive than the one in (1), unless one can 
identify a few states carrying a predominant contribution of the continuous 
part in (2). As shown in a simplified model [17], such a discretization would 
probably lead to the same numerical difficulties as above, as the necessary 
truncation of the integral in( 2) is not likely to be less expensive than a sum 
with large N max in (1). 

As a matter of fact one knows that the so-called metastable or resonant states 
corresponding to complex singularities of the resolvent give important con- 
tributions to the density of states. These states decay exponentially in time, 
preventing them to be considered as part of the spectrum, although they cor- 
respond to long-living nuclear configurations. Moreover it has been soon recog- 
nized (see [18] and references therein) that, despite their pathological asymp- 
totic behaviour, these unstable states could be used as generalized eigenfunc- 
tions that can be used to describe resonance phenomena in nuclear collisions. 
In this framework, a given state is expanded on a set of functions that includes 
resonant functions on the same footing as bound state eigenfunctions. 

The drawbacks of such a procedure are well known: as the resulting problem 
looses its self-adjointness, the complex spectrum has to be reinterpreted to- 
gether with the associated wave functions. However such expansions have been 
frequently used to describe various physical situations (potential and obsta- 
cle scattering, nuclear reactions... [19]). Recently, this approach has also been 
used to treat the problem of multiconfiguration mixing within the nuclear shell 
model [20]. 

Our purpose is to show that this old idea may be applied in a simple way 
to self-consistent computations when loosely bound or metastable states are 
expected, which is precisely the case in HFB computations for heavy nuclei, 
especially those which are close to drip-lines. 

This paper is organized as follows. In Sec. 2, we consider a model where the 
Schrodinger equation is solved in a bounded domain, and we apply it to a sim- 
plified one-body problem. In Sec. 3 some numerical experiments are presented, 
while Sec. 4 briefly presents our conclusions. 



2 The problem in a bounded domain 

Classically, in order to derive generalized eigenfunction expansions such as (2), 
one first computes the Green function associated to the problem and one uses 
complex analysis to continue this Green function and study its singularities. 
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Then, by applying the Cauchy residue theorem, one obtains (2), where the 
sum corresponds to the discrete spectrum, and the integral takes into account 
the continuous part of the spectrum (see Appendix A). 

A natural way to avoid the delicate discretization of the continuous contri- 
bution in (2) is to restrict the initial problem to a bounded region Q, where 
one expects that most of the physics takes place, and to impose "transparent" 
boundary conditions reflecting the properties of the exterior (unbounded) do- 
main. In fact, this last point deserves a special comment. It is clear that the 
kind of artificial conditions we impose on the boundary must be of particular 
nature in order to mimic the behaviour of the "exterior world", in a trans- 
parent way: the waves produced in the interior region must not be disturbed 
(reflected) by the boundary. For example, a standard Dirichlet or Von Neu- 
mann boundary conditions would lead to a problem essentially different from 
the one we want to solve and would produce results without clear connection 
with the solution of our original problem. One can check easily in particular 
that Dirichlet condition would lead to nothing but a Fourier series expansion, 
which makes the obtained solution strongly domain-dependent. 

One can show that this kind of "transparent" boundary conditions, which has 
been developped in the acoustic domain and more recently in the Schrodinger 
context (see [21] and references therein) leads to a complex condition pro- 
ducing a non-self adjoint problem: this is the price to pay for "replacing the 
exterior by the boundary". 

By restricting the domain, one can hope to get improved asymptotic properties 
for the modified Green functions, allowing one to bypass the integral along the 
physical cuts in (2). As a consequence, complex isolated singularities appear 
outside the imaginary axis, coming from the artificial complex conditions at 
the boundary of Q. 

To our knowledge, this idea goes back to Kapur and Peierls [22] who first intro- 
duced in the thirties non self-adjoint flux conditions associated to Schrodinger 
equations, a method which has been widely used in the context of nuclear 
reactions [18]. 



2.1 The model 

In the particular situation of the square well potential, the above method 
can be worked out in a very simple way. The idea is to replace the harmonic 
oscillator by a new reference potential V such that the hamiltonian H = T + V 
has a continuous spectrum, in order to describe more accurately the loosely 
bound states and, if necessary, metastable states. 
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As we expect that the physics takes place in the region where the potential is 
non-zero, we restrict the dynamics to fl and we impose a boundary condition 
at the boundary which mimics the external scattering behaviour. 

To be specific, let us consider the Schrodinger equation for a particle in a 
compactly supported spherically symmetric potential V: 

(--^ + V(r)J<p k (r) = k 2 <p k (r), for r G SI = (0, r ). (3) 

Where we have considered ft 2 /2m = 1 and an orbital momentum £ = 0. As a 
new reference potential V, we take the square well defined as: 



V(r) 



-V , for r < r , 

(4) 

0, otherwise, 



where V > and r > 0. 



We chose the square well because the associated eigenfunctions are rather 
elementary (Bessel functions). Moreover, in order to simplify the arguments 
we restrict the analysis to s-waves, in which case the eigenfunctions are simply 
sine and cosine functions. 

As an extra boundary condition for r = r , we choose the free radiation 
condition: 

(jL ~ ik ) = °' for r = r °' ^ 

From the above discussion about "transparent" boundary conditions, we stress 
that the non- self adjoint boundary condition (5) amounts to fix an imaginary 
flux condition reflecting exactly the physical radiation behaviour at infinity 
for the free Schrodinger equation. We feel that this choice corresponds to 
our purpose to produce the smallest perturbation near the artificial boundary 
r = r . 

The problem (3)-(5) is easily solved by computing the corresponding Green 
function (see Appendix A), giving simple explicit expressions for the associated 
eigenfunctions and the following expansion for any function / G L 2 (0,r ): 

f( r ') = ~~T — ( f sinp n rf(r) dr) sinp n r', (6) 
„ % + k n r \Jo J 

where the k n are the solutions both, real and complex, of the equation : 

p n cot p n r = ik n , (7) 
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with p n = (k^ + Vq) 1 / 2 . For real roots we recover bound and virtual state eigen- 
functions, while for non real ones we get resonant and anti-resonant eigenfunc- 
tions. 



2.2 Application to HFB problems 



Let us consider the familiar Hartree-Fock-Bogoliubov equations: 



f T + U V W^i 
v V -T-U) [ip 2 



E + X W^i 
E-XJ L 2 , 



where the particle-hole and particle-particle potentials U and V depend on ipi 
and ip2- 

In order to solve this non linear eigenvalue problem we expand the unknown 
wave functions ipi and ip2 on our reference basis, suitably truncated up to a 
rank N, depending on the desired accuracy: 



where 



N N 

^1 = a n¥n, and 1p 2 = briVn, 
n=0 n=0 



(p n (r) = sm p n r 



(9) 



(10) 



and the coefficients p n are obtained by solving the transcendental equation 
(7). 

Replacing the functions ipi(r) and ^(f) by there expansions, multiplying on 
the left by <p m (r) an d integrating, one gets the following expressions for (8): 



( N 



N 



^ ' (T mw -|- U mn )tt n -|- Vmn^n {.E -\- A) ^ ' R mn (l m , 
n=0 n=0 
N N 

"j V mn Q, n (T mn -|- U mra )6 n (E A) ^ ' R mn & TO . 



Kn=0 

or in matrix form: 
( 



71=0 



V 



T + U V 
V -T-U 



(E + A)R 
(E — A)R 



;i2) 



where T, U, V and R are complex hermitian matrices and \I/ is the column 
vector of components q = (ai, a^, &!, 6^). 
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The equation (12) is a generalized eigenvalue problem, from which the coeffi- 
cients of the expansion (10) together with the eigenvalues E can be obtained. 



If we solve the problem using only the functions <p n (corresponding to £ = 0), 
these elements are rather simple and in many cases they can be computed 
analytically. 

Some care must be taken in the evaluation of the kinetic energy T. As dis- 
cussed before, the expansion we use is valid only for r e (0, r ) and expansion 
(9) concerns only the part of the functions ^1,2 ( r ) defined in this interval. 
Consequently equations (9) have to be replaced by: 

Y(r - r)^i(r) = ^a n y? n (r), and Y(r - r)i/j 2 (r) = b n <p n (r), (13) 

n n 

where Y(r — r) is the Heaviside function. Taking the second derivative of this 
function, one gets: 



dr 2 



Y(r - r)ipi(r) 



dr 2 



dr 



r=r 



(14) 



As a consequence, the matrix elements of the kinetic energy T are: 



h 2 r r o 



<Pn(r) dr 



n -R mn + sin(/c m r ) cos(/c„r ) 



2m 



2m 



(15) 



They are computed in the same way as those of R (see relation (11)). Let 
us note that similar formulas have to be employed if the potentials U and V 
contain velocity-dependent contributions for example in the spin-orbit com- 
ponents. 

In the case of a nucleon-nucleon force of simple form (Skyrme or Gogny forces), 
the matrices U and V can be calculated analytically. Let us note that since T 
and R are hermitian matrices, the reality of the eigenvalues in equation (11) 
is insured. Finally the Coulomb potential can also be calculated analytically. 



2.3 Canonical states 



In order to interpret the HFB results, let us introduce the normal reduced 
density p(r, r')\ 
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rr'p(r,r') = ( r ) V4 ( r ')> (16) 

k 

where ip^ denotes the lower component of the k th solution of equation (12). 
With a similar notation, the abnormal reduced density p(r, r') is given by: 

rr'p(r,r') = -]T#(r) ip[ k) *(r'), (17) 
k 

Using expansion (9), these relations become: 

rr'p(r,r>) = ££&?°&f Vi(r) ^(/), (18) 

k i,j 

and: 

rr'p(r, r') = - b\ k) af ] ipi(r) <f*{r'), (19) 

k i,j 

for < r, r' < r . 

The interpretation of these quantities is clear (see for example [1]). The normal 
density is localized as soon as the chemical parameter A is negative. The 
canonical basis {ipj} is defined as the set of eigenvectors of the normal density: 

J dr'rr' p(r,r') ipiij-') =vftp i (r), (20) 

where the eigenvalues v f are the occupation probabilities of the corresponding 
canonical states. It is straightforward to show that if the are expanded in 
the same way as the quasi-particle wave functions, this equation corresponds 
to: 

pR*i = v^i, (21) 

where is the vector built with the coefficients of the expansion of ipi(r) and 
p is the matrix defined by: 

ft,=E?f ( 22 ) 

k 

It is important to stress that pij is not properly speaking a representation of 
the normal density since the set of functions {ipi} does not form a basis. 

Finally, we define the energies of the canonical states as the diagonal ma- 
trix element of the Hartree-Fock Hamiltonian in the basis of the canonical 
states [1]: 

e, = (%\H\%) (23) 
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This term can be easily written once the coefficients of the expansion of the 
canonical states are known. 



3 Numerical simulations 



To illustrate this method we have chosen to solve a simplified problem of HFB 
type with spherical symmetry. We consider the following equations: 





(24) 



where T is the kinetic energy operator, U and V are respectively the particle- 
hole and particle-particle mean-fields. We do not take into account the possible 
effects of a Coulomb field or of a spin-orbit potential, although this would not 
lead to any difficulties in the treatment. 



3.1 Hartree-Fock case (no pairing) 



We first consider the case without pairing, i.e. V — 0. In this case, the Lagrange 
multiplier A is not needed any more and we can omit it in the equations. We 
want to check the ability of our method to construct the correct bound and 
resonant states associated with a given potential. To this aim, we consider the 
gaussian potential used in [23]: 

W(r) = 5e- - 25 ^- 3 - 5 ) 2 -8e-°- 2r2 (25) 

This potential is given in MeV and the radial coordinate in fm in units where 
h 2 /2m = 1/2 MeV.fm 2 . This potential has a rich spectrum in all partial waves 
and the discretes states have been calculated with high accuracy in [23]. The 
energies and widths of the discrete states up to 10 MeV are reported in table 1. 

To build the set of functions (10) we have used a rectangular well of depth 30 
MeV and radius 12 fm. By solving equation (7) we have found 30 bound states 
and 29 virtual states. We complete this set with the lowest resonances and 
companion anti-resonances. The corresponding functions are not orthogonal 
and we cannot build the matrix R using all of them, otherwise the generalized 
eigenvalue problem (12) would be - at least numerically - singular. In order 
to reject some functions we use the following procedure: we build the full 
matrix R, using all functions and try to invert it with a simple Gauss method. 
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1 


E 


r 




e 


E 


r 





-4.571 183 







2 


-0.759 532 








-0.884 281 







2 


2.384152 


0.000 083 





2.252 381 


0.000118 




2 


4.659 669 


0.217375 





4.500 948 


0.247951 




2 


6.163 324 


2.408 615 





6.008 281 


2.516116 




3 


1.009 032 


~ 





7.587937 


6.266 307 




3 


3.810 922 


0.008196 


1 


-2.619 884 







3 


5.581406 


0.932 281 


1 


0.807635 


~ 




3 


7.091 723 


4.003469 


1 


3.577387 


0.013 017 




4 


2.676 525 


0.000 028 


1 


5.312 240 


1.063163 




4 


5.025176 


0.148 085 


1 


6.845 729 


4.224 511 




4 


6.522 260 


2.137256 


1 


8.427005 


8.441 884 




4 


8.059 380 


5.759133 



Table 1 

Discrete levels associated with the potential defined by (25). The energies and widths 
are determined by looking at the position of the poles of the S matrix in the complex 
energy plane. 

During this calculation, if we encounter a pivot smaller than a given value, 
then we reject the corresponding function. As a result, most of the times, we 
cannot use simultaneously bound and virtual states, because those two sets 
span spaces which are almost identical. In this example, we have been led 
to keep the 30 bound states, 3 resonances and 3 anti-resonances. The total 
number of functions used to expand the solutions is then 36. 

Because of the absence of pairing field, the matrix equation (12) is equivalent 
to the Hartree-Fock problem: 



T + U 



(26) 



The matrix elements of the matrices T, R and U can be computed analytically. 
Their expressions are given in Appendix B. After computing the matrices in 
(26) and diagonalizing the matrix Rr^T+U) we obtain a set of eigenvalues 
and eigenvectors that we normalize to 1 inside the domain (0,r ). 

We have reported in table 2 the eigenvalues obtained below 10 MeV. We can 
see that the values from table 1 are nicely reproduced. The bound states are 
given with an accuracy of five to six significant digits for the £ = and £ = 2 
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f — 1 


P — 9 


— Q 
-c — O 


P — A 


-4.571 182 


-2.605 981 


-0.759 533 


r\ 1 or* ono 

0.186 392 


0.229 617 


-0.884 280 


n 1 on on A 

0.130 394 


n 1 r - n r 

0.153 058 


n roc ,i n 

0.535 440 


n r r 1 on 

0.585 189 


Alio nnn 

0.118 909 


n a 71 on /i 

0.471 894 


n /i nT cr ot 

0.497 537 


1.009 141 


1 1 on 1 1 c 

1.120 116 


n a r non 

0.458 98U 


n <^ a a 

0.836 844 


1 noo ono 

1.023 303 


1.065 186 


1 onn nnn 

1.800 000 


n non rnn 

0.980 599 


n nn r 1 r T 

0.995 157 


1.692 437 


1 too nor' 

1.738 986 


r no 'i ao 

2.592 748 


1 c A A noo 

1.644 922 


1 f*f*c\ n a n 

1.660 949 


2.384 031 


r /i 00 /i 

2.524 324 


2.676 568 


2.252 302 


2.436 584 


O A TO 1 O /I 

2.472 184 


t n a n 

3.387 049 


O 1 /"* C A n CT 

3.466 495 


2.418 722 


noo 1 n 

3.283 319 


O OOT iVW 

3.327 296 


3.811 796 


A T /i r 

4.374 556 


otf cr 1 or* 

3.265 180 


3.617 422 


A 1 no om 

4.193 201 


A onn 010 

4.290 813 


/i n n non 

4.993 930 


4.106 893 


J 1 Tn /ion 

4.170 489 


4.675 124 


r* 1 on 

5.120 114 


r ono nn 1 

5.392 091 


4.546 091 


4.968 645 


r* 1 r onn 

5.215 890 


5.742 917 


6.207 120 


5.140 059 


bt r* n cr -1 /< -1 

5.605 141 


6.024 047 


c acts, r nn 

6.469 599 


i-f nno oot 

7.002 387 


5 946 712 


6 379 345 


6 841 667 


7 342 215 


7 878 341 


6.773 633 


7.256 568 


7.743 484 


8.286 709 


8.849 874 


7.684636 


8.230 319 


8.731860 


9.300 691 


9.897122 


8.679 769 


9.252138 


9.793 386 






9.746 902 











Table 2 

Eigenvalues obtained by diagonalizing the matrix equation (26). The boldface num- 
bers correspond to the energies of the bound states and of the scattering states 
strongly localized inside the potential. 

states. The narrow resonances are also very well reproduced. The discrepancies 
between the values in table 1 and 2 for the resonances energies are mainly due 
to their widths. The wide resonances found in table 1 by computing the S 
matrix cannot be seen directly in table 2 because they spread over a number 
of scattering states. The states with £ = 1 are not so well reproduced, but 
we have checked that the agreement is better if we use more functions to 
expand the solutions (in this example, all the results are obtained by using 
36 functions to expand the solutions). A plausible explanation concerning this 
deficiency is the following: the set of functions used to expand the solutions 
behaves like r in the limit r — > 0, so the convergence is supposed to be slower 
for partial waves with £ > 0, the most important consequences are seen on 
waves with I = 1 because these functions have a significant amplitude near the 
origin where they are not supposed to be correctly reproduced. This effect is 
smoothed for higher partial waves because the centrifugal term prevent them 
to "explore" the vincinity of the origin. It is really remarkable to see that the 
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1 


A E = -4.571 182 MeV 

\ 




E = -0.884280 MeV 


- 

-0.5 











2 4 6 R 10 12 

r[fm] 


2 4 6 8 10 
r [fm] 




E = 2.252302 MeV 




E = 4.546091 MeV 













2 4 6 8 10 12 

r[fm] 


2 4 6 8 10 
r [fm] 



Fig. 1. Wave functions obtained for £ = 0. Upper-left and Upper-right: the two 
bound states at -4.571182 MeV and -0.884280 MeV. Lower-left: the first resonance 
at 2.252302 MeV. Lower-right: the scattering state at 4.546091 MeV, this state is 
not much localized but can be identified as a resonance. 

solutions with orbital angular momentum from to 4 can be reproduced by 
expansions on functions with £ = 0. 




Fig. 2. Wave functions obtained for i = 1. Right: wave function of the bound state 
at -2.605981 MeV. Middle: the first resonance at 0.836844 MeV. Right: second 
resonance at 3.617422 MeV. Small unphysical oscillations can be seen on the wave 
functions, this phenomene is related with the behavior of the I = 1 functions near 
the origin, see the text for more detailed explanations. 

We can also observe how the eigenfunctions behave. Some of them with £ = 
0, 1 and 2 are plotted on figures 1, 2 and 3. Figure 1 shows that the method is 
able to reproduce the deeply bound states, weakly bound states and scattering 
states. In the particular case of the scattering state at 2.252302 MeV, the wave 
function is strongly localized in the region of the potential and corresponds to 
a resonant state. 

On figure 2 we see that the wave functions for £ = 1 are less correctly repro- 
duced. Some small unphysical oscillations appear on the tail of the functions. 
As it was discussed previously, this problem can be cured by enlarging the set 
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E = -0.759533 MeV 



E = 2.384031 MeV 





02 468 10 12 02468 10 12 0246 8 10 12 

r [fm] r [frn] r [fm] 

Fig. 3. Wave functions obtained for I = 2. Left: Bound state at -0.759533 MeV. Mid- 
dle: Resonant scattering state at 2.384031 MeV. Right: Resonant scattering state 
at 4.675124 MeV. The second resonance is at higher energy so the wave function is 
less localized in the interior of the potential than the first one. 

of functions used to expand the solutions. 

We can see on figure 3 that despite the fact that the expansion of the solutions 
involves only functions with £ = 0, the behavior of the £ = states at the 
origin ~ r i+1 ) is very well reproduced. One observes that the 2.384031 
MeV resonance wave function is strongly localized in the region where the 
potential is attractive, while the 4.675124 MeV one is much less localized. 



3.2 Hartree-Fock-Bogoliubov case 



We choose U to be a Saxon- Woods potential: 

W 



W(r) = 



1 + e~ 



R 



(27) 



As the pairing field is known to be mainly localized around the nuclear surface, 
we adopt for the particle-particle channel potential V the derivative of a Saxon- 
Woods potential with the same radius and diffuseness as in the particle-hole 
channel, but with a different intensity: 



V(r) 



d_ 

dr 



Vn 



-Rq 



.1 + e~ 

The system (24) has been solved for partial waves from £ = to 4. 



(28) 



As discussed in [3], the spectrum of (24) is unbound from above and from 
below. The solutions for negative and positive energies are related by: 



M'E,r) = -ip 2 (E,r) 



(29) 



Therefore, when solving equation (24), we consider only the solutions with 
E > 0. 
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The numerical parameters of the particle-hole potential (28) have been taken 
as 



U = 32 MeV, r = 3.7 fm, a = 0.65 fm (30) 

and those of the particle-particle channel potential: 

V = 4 MeV, r = 3.7 fm, a = 0.65 fm (31) 

The energy scale is chosen such that ft 2 /2m = 20 MeV in this example. 

The set of functions {ipi} used to expand the solutions was generated from a 
box of radius 40 fm and depth 180 MeV which contains 38 bound states and 
37 virtual states. The set is built with the 38 bound states and the 2 first 
couples of resonances and anti-resonances, so the dimension of the matrices is 
42. 



e 








1 


2 


E (MeV) 


-19.288 


-0.856 


-9.460 


-0.124 



Table 3 

Energies of the bound states in the Saxon- Woods potential (without pairing). 

In table 3 we have reported the energies and the norms of the particle-particle 
bound states when the pairing field and the chemical parameter are set to 
zero. The system has one loosely bound state 124 keV below the continuum. 
As mentioned in the previous paragraph, the convergence in the basis of the 
different partial waves is quite fast and good except for £ = 1 for which it is 
slower. For this reason, we will present in detail the results only for £=0 and 
2. As in the previous paragraph, we have checked that a satisfactory conver- 
gence for £ = 1 can be obtained by enlarging the set of functions {</?«}. In a 
future development of this method, where most general case of 3 dimension- 
nal computation will be investigated, this problem will not appear because 
the effective centrifugal potential is not present anymore. 




E = 0.936 MeV 



E = 18.594 MeV 



Fig. 4. Wave functions of three quasi-particle states for I = 0. The two components 
ipi (solid line) and tp2 (dashed line) have the same asymptotic behaviour for the 
state on the left panel which corresponds to a discret level (E < —A), while different 
asymptotic behaviours are obtained for the two others which correspond to states 
with quasi-particle energies E > — A. 
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In table 4 we have reported the energies and the norms N 2 of the lower com- 
ponents of the quasi-particle states for which N 2 is larger than 0.001, for £ = 
and £ = 2. For higher orbital momenta, the lower component of all solutions 
is almost zero and does not contribute to the particle density. With the value 
chosen for the chemical potential (A = —0.750 MeV), the number of parti- 
cles in the system is on the average 3.163 for £ = and 0.417 for £ = 2. We 
have represented in figure 4 three quasi-particle wave functions: two states 
close to the Fermi sea and a deep hole state, with quasi-particle energies 0.475 
MeV, 0.936 MeV and 18.594 MeV, respectively. The energy of the first quasi- 
particle state is smaller than |A| and therefore this state belongs to the discrete 
spectrum (this is actually the only discrete state in the full spectrum). The 
asymptotic behaviour of the quasi-particle wave function is [3]: 



ipi(r) oc exp(— Kir) 
^ 2 {r) oc exp(-K 2 r) 



for r — > 00. 



(32) 



The two other quasi-particle states shown in the figure are continuum states 
with an asymptotic behaviour given by: 



ipi(r) oc sin(A; 1 r) 
^2(7") oc exp(—K 2 r) 
with obvious definitions for k±, K\ and k 2 . 



for r — > 00, 



(33) 



We can see in figure 4 that the method we employ here nicely reproduces those 
two completely different asymptotic behaviours. Using the relation given in 





I ~- 


= 




£ = 


2 


E (MeV) 


N 2 


E (MeV) 


N 2 


E (MeV) 


N 2 


0.475 


0.549 


4.662 


0.001 


1.072 


0.152 


0.936 


0.013 


15.338 


0.001 


1.157 


0.043 


1.443 


0.010 


18.212 


0.090 


1.742 


0.004 


2.239 


0.005 


18.594 


0.905 


2.587 


0.002 


3.315 


0.002 


21.402 


0.001 


3.689 


0.001 



Table 4 

Energies E and norms N 2 of the lower components of the quasi-particle states with 
orbital momenta £=0 and 2. Only the states with N 2 > 0.001 have been reported. 

equation (16) the densities of particles for £=0 and 2 can be built. These 
quantities are represented in figure 5 using linear and logarithmic scales. 

We notice on the logarithmic plots that the behaviour of the density becomes 
oscillatory when it falls below 10 -9 fm -3 . This is due to the truncation of the 
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basis expansion and also to some numerical inaccuracies. However the particle 
density is so small in this region that this pathology should not play any 
significant role in the calculation of any observable. 




r [fm] r [fm] 



Fig. 5. Particle densities in fm 3 for t = (top) and i = 2 (bottom) as functions 
of r. Figures on the right hand side are in logarithmic scale. 

Eventually we diagonalize the density operator in order to construct the canon- 
ical states. The single-particle energies and occupations of £ = canonical 
states are displayed in table 5. The wave- functions of the first three ones are 
shown in figure 6. We notice that the expected decrease to zero at large r of 
these wave functions is remarkably well described within the present approach. 



e (MeV) 


-18.570 


-0.127 


5.919 


24.570 


v 2 


0.9997 


0.5802 


0.0012 


0.0002 



Table 5 

Energies e and occupations v 2 of the canonical states obtained for I = 0. Only the 
states with an occupation greater than 0.0001 have been reported. 




r [fm] r [fm] r [fm] 



Fig. 6. Wave functions and energies of the three first canonical states for I = 0. 
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4 Conclusion 



We have presented in this work an expansion method which can be used to 
express both discrete and continuum solutions of the most general HF or HFB 
equations, when any kind of nucleon-nucleon interaction is used. Because of 
the simplicity of the functions on which we expand the solutions, many steps 
of the resolution can be done analytically This is the great advantage of the 
method. For example, the matrices needed for the resolution of the equations 
(12) can be evaluated analytically if the mean-field potentials are of Saxon- 
Woods, and derivative of Saxon- Woods or Gaussian forms. 

The fact that this method avoids the problems associated with the discretiza- 
tion of integro-differential equations has an important consequence: it is fast 
and accurate. The only part of the problem that has to be done numerically is 
the inversion of the matrix R and the diagonalisation of (12). In general, the 
dimension involved are rather small, at least in the case of spherical symme- 
try. So these steps can be done rapidly and with a good accuracy. In addition, 
as the quasi-particle states are expressed as finite expansions on a basis of 
wave-functions, the method ensures that the variational HF and HFB proce- 
dures lead to an upper bound of the total energy of the system, which can be 
improved at will by increasing the size of the basis. 

Because of the form of the kinetic energy (15), the matrix problem that we 
solve is formally equivalent to solving the corresponding integro-differential 
equations in a box of radius tq with the non local boundary condition \I/'(ro) = 
iK^/(r Q ) for the solutions, where K is the operator defined by: 



and k n are the solutions of equation (7). In this sense, the present method is 
completely equivalent to working directly in the coordinate representation as 
done e.g. in [3] and, provided the size of the box is the same, comparable results 
would be obtained. The advantage of the technique we propose is that it can 
be applied to HFB calculations employing a finite-range nuclear interaction 
such as the Gogny force. Actually, once the matrix elements of the two-body 
force are computed, the method allows one to solve the HF and HFB equations 
as easily as with a basis of harmonic oscillator states, which has been done 
extensively with the Gogny force [2,14]. In addition, extensions of the method 
to non-spherical nuclei could be done in the same fashion as with harmonic 
oscillator bases. 

Let us note that in the limit where the depth of the rectangular well goes 
to infinity (Vb — > +oo in equation (4)) usual discrete Fourier expansions are 
recovered. In that case the formalism is greatly simplified because the tran- 
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scendental equation (7) becomes trivial and the overlap matrix R reduces to 
unity. 



The more general case investigated here has however the advantage of defining 
a reference set of functions which are adapted to the expansion of the (quasi- 
)particle state solutions of the HFB equations in nuclear systems. As shown 
in the present work, good approximation of the solutions can therefore be 
obtained, using comparatively small expansions. 

The implementation of the present technique in the fully self-consistent HFB 
procedure will be the next step of this study. An important issue will be to 
check if nuclear properties can be accurately reproduced in realistic situations, 
in particular in nuclei close to drip lines, by using bases small enough to ensure 
extensive calculations even in heavy nuclei. If such a requirement is met, the 
method of solving the HFB equations we propose should open a broad range 
of new nuclear studies for the future. 



A The decomposition of identity 



In this paragraph, we consider ti 2 /2m = 1 to simplify the expressions. We con- 
sider the Green function G KP (r,r';k) (KP for Kapur-Peierls), corresponding 
to mixed boundary conditions on the domain ]0,ro[, solution of the (non self 
adjoint) problem: 



< 




(0,r'; fc) = 0, 



+ (V(r) - k 2 ) G KP (r, r'; k) = 5(r - r'), if < r, r' < r , 



(A.l) 



d 

^ dr 



G KP (r,r'; k)\ 



ikG KP (r, r'-k)\ 



If we put p = \fk 2 + V , the solution reads: 




tpr , for < r < r', 



(A.2) 



G^(r, r'\ k) = A 2 e ipr + B 2 e~ ipr , for r' < r < r , 
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The coefficients are given by the boundary conditions: 

GW(0,r';k) = 0, 
G^(r',r';k) = G ( ^(r',r';k), 
(G^)'(r' + , r>; k) = (G^)'(r'_, r'; k) - 1, 
(G (2) )'(ro, r'; fc) = ikG^(r , r'; k), 



(A.3) 



and we obtain: 



G^ p (r,r';£;) 



ifcsmp(r — r ) +pcosp(r — r ) smpr . 

1 , it < r < r , 

ik sinpr — pcospr p 

ik sinpfr — r ) + p cos p(r — r ) sinpr' . 

; , it r < r < r , 

ik sin pr — p cos pr p 



(A.4) 



One checks that k — > G KP (r,r'; k) can be meromorphically continued in the 
complex plane and that for < r,r' < r , it decays exponentially when 
\k\ — > oo in the whole complex plane C. 

In order to get a completeness formula, we consider a compactly supported 
C 2 function /, (which is zero near the end points r = and r = r , and we 
denote by g the function defined by: 



9(r) 



-Vo)f(r). 



(A.5) 



By multiplying (A.l) by /, (A.5) by G KP (r, r'\ k), subtracting the two resulting 
equations and integrating over (0,r ), we get: 

kG KP (r, r'; k)f(r) dr = -±f{r>) + \ jT G KP (r, r>; k)g{r) dr. (A.6) 

By integrating on each side on a big circle Ca = {k e C : \k\ — A}, we 
obtain: 



lim — / [ f r ° kG KP {r, r'; fc)/(r) 

fe-*oo 2«7T JC A Uo 



dk = -f(r>). 



(A.7) 



By applying the Cauchy residues theorem, we get the following formula: 



r pro 

f(r') = -Y^R es / kG KP (r, r'; k)f(r) dr 
„ Uo 



(A.8) 
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where the k n are the solutions of the equation : 

p n cot p n r = ik n , (A. 9) 



with p n = yk^ + Vo After computing the residue, we get the following decom- 
position identity: 

f( r ') = J2 —T — / sin Pn.rf{r) air sinp n r'. (A.IO) 

In this identity, the index n labels all the complex solutions of (A. 9). 

One may assume for convenience that the first A" roots are pure imaginary 
k n = iK n with K n > 0, and correspond to the A" real eigenvalues E n = —n n , 
1 < n < N, with —Vq < E n < 0. The other roots are located symmetrically 
in the negative half-plane ^sm{k) < 0, in such a way that if k is a root, —k is 
also a root of (A. 9). 

Numerically, the pure imaginary roots k n can be computed by using Newton's 
method. For the other complex roots, one can show that [24] the behaviour of 
in = k n r = x n + iy n for n large, is given by: 

1 , f 2c n 

Xn = C n lo §(^- + -» 

° n (A.11) 
^ = -log(^)-^(log(^)-l) 2 + ... 

where c n = (n + |)7r. Then £ n ~ nix — i log n gives an initial guess for the large 
order roots. 



B Matrix elements of the potentials: 



Gaussian Potential 



To express the matrix elements of the Gaussian potential in compact form, 
we introduce the function: 



2 

exp(/37 + j^K/jr 
8^a 



erf 



v^-/3-^)J+erf(v^(/3 + ^ 



(B.l) 
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where erf denotes the error function [25]. The matrix elements of a general 
gaussian function: 

-a(r- (3) 



fae(r) = e-«r-f>, (B.2) 



can be written: 

-L 



J sm(K*r)f afS (r) sin^r) dr = [i(k* + re,)) + + /%) 

- - «,-)) - - «*)). 

(B.3) 

If we choose (ct = 0.25, /3 = 3.5) and (a = 0.2, /3 = 0), we can use this relation 
to compute the matrix element of the potential (25) very simply. 

Saxon- Woods potential 



In order to compute the matrix elements of the Saxon- Woods potential, we 
introduce the function i^ (/c) defined by: 



) = / T^dr. (B.4) 

J0 1 + e a 



I ro J K 



This integral can be evaluated using the hypergeometric function 2 F\ [25]: 



2F1 (l, iaK; 1 + iaK; — e ^ 

B.5) 

e iftL 2 Fi (l, ia/c; 1 + iaK; -e^) 



The matrix element of a Saxon- Woods potential V(r) with depth Vo, diffuse- 
ness a and radius r are given by: 



V t/ /- L sin«r) sin(K w r) 

V mn = - Vo ^ or, (B.6) 

JO X _|_ g a 

and can be written explicitly using the /£ a (/c) functions: 
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2 



Lr a\ K m K n) 



Ir a( K rn + K «) 



(B.7) 



If K* m = Kn, this expression becomes: 



2 



L + a In 



1 + ev 

_L Itt 

e« + e a 



(B.8) 



Derivative of Saxon- Woods potential 



We have used for the particle-particle channel a potential with the shape of 
the derivative of a Saxon- Woods potential with depth A , diffuseness a and 
radius r : 



A(r) = A 



d 



dr 



1 + e 



(B.9) 



The matrix elements of this potential can also be written using the func- 
tion (B.5): 



Ao 



K m K n 



Ai 



1 r a\ K m ~ K n) — -^roaK^n ~ K r , 



+ 



sin(K,* n L) sm(K n L) n* m + K n 



l + exp(^) 



4i 



(B.10) 



Centrifugal term 



The matrix elements associated with the centrifugal term are given by: 



C 



h 2 £(£ + 1) f L sin /t*r sin Kf 



2m 



dr. 



(B.11) 
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This integral can be evaluated using the sine integral function [25]: 
h 2 £(£+l) 



Cij 



2m 



sin L sin k^L Ka Ka „. , * . r 
—= — h „ Si(k* - Kj)L 



^±^Si« + k,)L 



(B.12) 
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